################################################################################
########################     Analysis and Plots     ############################
################################################################################

library(Matching)
library(tidyverse)

################################################################################

rm(list=ls())
load("~/data/virginclip_pop.RData")

################################################################################

# 1. No trimming on forest cover in 2005

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW>0) 

# outcome vars

virginclip_pop$did2016 <- virginclip_pop$forestha_2016 - virginclip_pop$forestha_2011
#virginclip_pop$did2005 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2005

virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011
#virginclip_pop$did2004 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$did2018 <- virginclip_pop$forestha_2018 - virginclip_pop$forestha_2011
#virginclip_pop$did2003 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2003

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017
t4 = 2018

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2016 <- virginclip_pop$did2016 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2018 <- virginclip_pop$did2018 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                           forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                           distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                           agc_bcc_mean + pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

virginclip_pop$pscore<- ps_virginclip_pop$fitted

##### DDD REGRESSIONS #####

# 1:2 matching

PSM_DDD_virginclip_pop2017_12 <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                     X = virginclip_pop$pscore, estimand = "ATT", M = 2, ties = T,
                                     replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2017_12)


# 1:3 matching 

PSM_DDD_virginclip_pop2017_13 <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                     X = virginclip_pop$pscore, estimand = "ATT", M = 3, ties = T,
                                     replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2017_13)

# 1:5 matching

PSM_DDD_virginclip_pop2017_15 <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                     X = virginclip_pop$pscore, estimand = "ATT", M = 5, ties = T,
                                     replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2017_15)


# Elevation > 1000 m
rm(virginclip_pop, ps_virginclip_pop)
load("~/data/virginclip_pop_noelev.RData")

################################################################################

# 1. No trimming on forest cover in 2005

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW>0) %>%
  filter(elev > 1000)

# outcome vars

virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017
t4 = 2018

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                           forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                           distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                           agc_bcc_mean + pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

virginclip_pop$pscore<- ps_virginclip_pop$fitted


###### DDD REGRESSIONS ######

PSM_DDD_virginclip_pop2017_1000 <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                     X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                     replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2017_1000)


# All Elevation 
rm(virginclip_pop, ps_virginclip_pop)
load("~/data/virginclip_pop_noelev.RData")

################################################################################

# 1. No trimming on forest cover in 2005

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW>0) %>%
  filter(elev > 1000)

# outcome vars

virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017
t4 = 2018

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                           forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                           distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                           agc_bcc_mean + pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

virginclip_pop$pscore<- ps_virginclip_pop$fitted


###### DDD REGRESSIONS ######

PSM_DDD_virginclip_pop2017_ae <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                     X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                     replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2017_ae)


rm(virginclip_pop, ps_virginclip_pop)
save.image("~/results/table-s18.RData")



